Density inhomogeneities in heavy ion collisions around the critical point 
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We study the hydrodynamical expansion of a hot and baryon-dense quark fluid coupled to clas- 
sical real-time evolution of the long wavelength modes of the chiral field. Significant density in- 
homogeneities develop dynamically when the transition to the symmetry-broken state occurs. We 
find that the amplitude of the density inhomogeneities is larger for expansion trajectories crossing 
the line of first-order transitions than for crossovers, which could provide some information on the 
location of a critical point. A few possible experimental signatures for inhomogeneous decoupling 
surfaces are mentioned briefly. 
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I. INTRODUCTION 

Heavy ion collisions at high energies produce hot and 
baryon-dense strongly interacting matter and so pro- 
vide the opportunity to explore the phase diagram of 
QCD yj. Recent lattice QCD calculations at finite 
baryon-chemical potential 2] indicate that at sufficiently 
large baryon density a line of first order transitions ex- 
ists in the plane of temperature T versus baryon-chemical 
potential [is- This line separates the region where chiral 
symmetry is broken (as in vacuum) from that where it is 
approximately restored. Moving counter-clockwise along 
this phase boundary, i.e. towards higher T and lower /is, 
results in weaker first-order transitions and finally the 
line of first order transition ends at a second-order critical 
point [3J. Simulations with semi-realistic quark masses 
locate the endpoint at T E ~ 160 MeV, (i B ,E ^ 360 MeV. 
For /iB < I^b,e no phase transition in the strict sense oc- 
curs. Rather, the low- and high-temperature phases are 
continuously connected by a rapid crossover. 

Our goal here is to analyze the homogeneity of the 
"fluid" of QCD matter as it expands and cools. In par- 
ticular, we shall study expansion trajectories passing on 
either side of the critical point (i.e. either crossover or 
first order phase transition). As we will show, in the 
vicinity of the critical point the expanding fluid devel- 
ops significant inhomogeneities. Such density perturba- 
tions should also be present on the decoupling surface of 
hadrons. This is normally neglected in hydrodynamical 
simulations of heavy-ion collisions, which commonly as- 
sume that the hadrons freeze out at a fixed temperature 
or density. 

The fact that decoupling surfaces are typically not ho- 
mogeneous is well known. Take, for example, the WMAP 
data on the temperature fluctuations of the cosmic mi- 
crowave background (CMB) 0. The background pho- 
tons exhibit temperature fluctuations on the order of 
AT/T ~ 10~ 5 . From the CMB multipoles one hopes to 
gather information on their primordial origin. In heavy- 
ion collisions, on the other hand, we might get hints 
about the QCD phase transition, if it occurs shortly be- 
fore decoupling. The basic idea is that if a phase transi- 
tion occurs then it might leave imprints on the (energy-) 
density distribution on the freeze-out hypersurface. In 



particular, we expect the inhomogeneities to be smaller 
for crossovers and stronger for first order transitions. 

The first source of spatial density inhomogeneities in 
heavy-ion collisions is due to fluctuations in the number 
of participants and the number of collisions among the 
beam nucleons. Those fluctuations lead to an inhomoge- 
neous deposition of energy and of baryon number at cen- 
tral rapidity. From the event generator UrQMD, energy 
density inhomogeneities on the order of Ae ~ 1 GeV/fm 3 
have been predicted for central Pb+Pb collisions at top 
SPS energy (y/s ~ 17^4 GeV) which originate from fluc- 
tuations in "soft" (low momentum transfer) interactions; 
see fig. 4 in [|. At RHIC energies (y/s ~ 100-200A GeV) 
Ae could increase by up to an order of magnitude due 
to the additional semi-hard ("minijet") component, as 
discussed in ref. These authors also discussed the 
evolution of such initial-state inhomogeneities using equi- 
librium hydrodynamics in the ideal fluid approximation. 
They find that even the huge initial perturbations pre- 
dicted by the minijet model are strongly washed out un- 
til freeze out by the hydrodynamic expansion of the hot 
matter. Qualitatively, this can be understood from the 
following simple argument. The radius L ~ 1 fm of a 
hotspot grows linearly in time while its density drops in- 
versely proportional to its volume so that Ap ~ 1 /t 3 for 
expansion in three dimensions. The overall duration of 
the hydrodynamic expansion is expected to be (at least) 
on the order of the radius of the colliding nuclei p], Q , 
roughly ~ 5L. Therefore, any initial density concentra- 
tion should be diluted by about a factor of 100. Hence, 
initial perturbations of order one would leave traces on 
the percent level only, at the time of decoupling. For 
other studies of early-stage density inhomogeneities and 
their hydrodynamic evolution see refs. and , respec- 
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Here, we discuss a different source of density inhomo- 
geneities, namely those possibly generated in the course 
of a non-equilibrium transition from a (nearly) chirally 
symmetric state at high temperature and density to the 
state of broken symmetry at decoupling. The rather 
rapid transition expected to occur in high-energy heavy- 
ion collisions very likely forces the long wavelength modes 
of the chiral condensate out of equilibrium. This should 
then reflect in a rather non-uniform distribution of en- 
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ergy and baryon density in space (or, more precisely, on 
the decoupling hypersurface) . Most notably, since "freeze 
out" (decoupling of all particles) in heavy-ion collisions 
occurs shortly after the transition to the broken phase, 
such perturbations generated in the late stages of the evo- 
lution could largely survive and leave detectable traces in 
the final state. In this regard, we also recall the results of 
ref. 01 wn0 employed the collisionless Vlasov equation 
to study the real-time evolution of small initial density 
fluctuations within the NJL model. They observed an in- 
crease of the fluctuations already for the case where the 
expectation value of the chiral condensate was fixed to its 
equilibrium value, cf. their eq. (2). Here, we also treat 
the non-equilibrium dynamics of the chiral condensate, 
cf. our cq. Q below, which probes the structure of the 
effective potential in the vicinity of the critical point. 

II. THE MODEL 

For the current studies we extend the model from 
ref. |12| to allow for nonvanishing baryon density p. The 
Gell-Mann-Levy Lagrangian ^3] 

£ = q [«7 M <9 M - g(a + 75T ■ 7?)] q 

+ \ {d»v? + \{d^) 2 -U{a^) . (1) 

provides an effective theory for chiral symmetry breaking 
in QCD. It describes the interaction of two flavors of 
constituent quarks q = (u, d) with the chiral field <j) a — 
(a, 7?) . The potential, which exhibits both spontaneously 
and explicitly broken chiral symmetry, is 

U(a,7r) = ^{a 2 +Tf 2 -v 2 ) 2 -h q a-U . (2) 

The vacuum expectation values of the condensates are 
(cr) = /„. and (7?) =0, where f n — 93 MeV is the pion 
decay constant. The explicit symmetry breaking term 
is due to the non-zero pion mass, h q — f^rrQ, where 
tti-jy = 138 MeV. This leads to v 2 = f 2 - m\j\ 2 . The 
value of A 2 = 20 leads to a er-mass, m 2 = 2X 2 f 2 + m 2 , 
approximately equal to 600 MeV. 

We assume that the quarks constitute a thermal- 
ized fluid, which provides an expanding background in 
which the long-wavelength modes of the chiral conden- 
sate evolve. Integrating out the quarks generates an ef- 
fective potential for (j) a ; computing to one loop and for 
a homogeneous background (on the scale of a "fluid ele- 
ment"), this contribution is given by 

V eS (<t>, T, m) = U(cf>) -d q T J -0^ (3) 
{log(l + e^/ T )+0^- M )} . 

Here, d q = 12 denotes the color-spin-isospin degeneracy 
of the quarks and fi = the quark-chemical poten- 

tial. The two terms inside the integral correspond to 



the thermal contributions of quarks and anti-quarks, re- 
spectively; a (divergent) vacuum contribution has been 
absorbed into the T and fx independent potential U. V c s 
depends on the order parameter field through the effec- 
tive mass of the quarks, m 2 — g 2 (f> 2 , which enters the 

expression for the single-particle energy E — \jv 2 + m q - 

For sufficiently small quark-chemical potential /i one 
finds a smooth transition to approximately massless 
quarks at high T . For larger chemical potential, however, 
the effective potential exhibits a first-order phase tran- 
sition 14] . Along the line of first-order transitions the 
effective potential exhibits two degenerate minima which 
are seperated by a "nucleation barrier" . This barrier de- 
creases with fi and the two minima approach each other. 
At finally, the barrier vanishes, and so does the la- 
tent heat. For g = 3.3, which leads to a constituent quark 
mass in vacuum of « 307 MeV, the second-order critical 
point is located at T E ~ 100 MeV, fi E ~ 200 MeV. In- 
creasing the quark-field coupling g moves the endpoint 
E towards the temperature axis |15| (/j,e becomes =0 at 
about g « 3.7 and to slightly higher temperature. 
In what follows, we fix g = 3.3. 

The location of the endpoint does not agree quanti- 
tatively with that from recent lattice QCD studies with 
realistic quark masses, which find Tg « 160 MeV and 
He ~ 120 MeV 0. This failure is common to sev- 
eral models, cf. fig. 6 in and could be due to the 
neglect of heavier resonance states in the above effec- 
tive Lagrangian, see e.g. the discussion in 0|. Also, 
if deconfinement and chiral symmetry restoration occur 
simultaneously then the energy density contributed by 
Polyakov loops in the deconfined phase should be in- 
cluded as well [is|. Finally, note that the model fails 
to describe the nuclear matter ground state, which has 
non-zero pressure and is located in the phase coexistence 
region (Fig.^below). Nevertheless, one might hope that, 
qualitatively, the dynamics of relativistic quark fluids 
near the endpoint is not affected by the deficiencies of 
this simple model. 

The classical equations of motion for the chiral fields 
are 

Wa + ^ = 0. (4) 

0(f) a 

We do not account explicitly for a damping term due to 
decay processes or elastic collisions of the particles form- 
ing the condensate |l9j |. An ensemble average over ran- 
dom initial field configurations implicitly introduces such 
effects at the classical level. The (in-)accuracy of this 
approximation should be a matter of further study |20| . 
In J3J , the only explicit damping of field oscillations arises 
from the expansion of the fireball. 

The dynamical evolution of the thermalized degrees of 
freedom (fluid of quarks) is determined by the conserva- 
tion laws for energy, momentum and (net) baryon charge: 

9, CCid + T) = 
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^ (pa") = . (5) 

Here, u^ 1 is the fluid four- velocity and T^ id its energy- 
momentum tensor, which we assume to be of perfect fluid 
form. T^ u , in turn, is the energy-momentum tensor of 
the classical fields which can be obtained from the above 
Lagrangian in the standard fashion 0, [2l| . Note that 
we do not assume that the chiral fields are equilibrated 
with the heat bath of quarks. Hence, the fluid pressure 
depends not only on the energy and baryon density in the 
local rest frame but also on the chiral (order-parameter) 
field, i.e. p = p(e, p, <fi). 

We employ eq. (@J to also propagate initial field fluc- 
tuations through the transition; that is, our initial con- 
dition includes some generic "primordial" spectrum of 
fluctuations (see below) which then evolve in the effec- 
tive potential generated by the matter fields. 

III. RESULTS 
A. Initial Conditions 



potential near the transition. Of course, larger fluctua- 
tions would amplify the effects shown below. We corre- 
late the initial field fluctuations over approximately 1 fm 
as described in Our focus is on how those "primor- 
dial" fluctuations evolve through the various transitions. 

For definiteness, we shall consider two different sets of 
initial conditions: for set (I) we start the evolution at 
fixed initial energy density e cq = 2.8eo but vary the ini- 
tial baryon density p cq = (0, 0.6, 1.6, 2.1, 2.4, 2.8)po; for 
set (II), on the other hand, we start at fixed initial net 
baryon density p cq — 1.7 po but vary the initial energy 
density e eq = (1.4, 1.9, 2.9)eo- Here, eo and po denote 
nuclear matter ground state energy and baryon density, 
respectively. For low baryon density (I) (high energy 
density (II)), the expansion will then proceed through 
a crossover, while a baryon dense (I) (energy dilute (II)) 
droplet will decayvia a first-order phase transition. (In 
contrast, in 0, the type and strength of the tran- 
sition was controlled via the coupling constant g rather 
than the baryon density.) Our goal is to analyze the 
evolution of baryon density inhomogeneities. 



We employ the following set of simple initial conditions 
to illustrate qualitative effects. At t — 0, we initialize a 
sphere of hot and dense quarks with radius R = 5 f m and 
no initial collective motion, v(t = 0) = 0. The energy and 
baryon-density distribution is taken as 



e(t = 0,x) 
p(t = 0,x) 



l + exp(^) 



l + exp(^) 



(6) 



with a surface thickness of a = 0.3 fm. 

Within that sphere, the average chiral field corre- 
sponds to the minimum of V e s (e oq , p cq ). Specifically, we 
choose 



a(t = 0,x) = 6a(x) + U + - 2S_ 

1 + exp 

= 6a(r,<p,6) + (a)(r) 

Tr(t = 0, x) — Sir , 



'r-R\ 



(7) 



with 



the expectation value of the a field cor- 



responding to e cq and p cq . Thus, the chiral condensate 
nearly vanishes at the center, where the energy density 
of the quarks is large, and then quickly interpolates to 
j v where the matter density is low. The system sub- 
sequently expands hydrodynamically on account of the 
nonzero pressure. 

Sa(x) represents Gaussian random fluctuations of the 
fields which are distributed according to 



P[ty«]oc exp {-Hl/2 {8(g)) 



(8) 



The results presented here were obtained with a width of 
y/(Sa 2 ) = v/3, y/(5if 2 ) = 0. These relatively moderate 
amplitudes suffice to probe the structure of the effective 



B. Time evolution 

Fig. shows the trajectory of the system within the 
phase diagram for both sets of initial conditions. For 
simplicity, we chose a foliation of space-time by flat hy- 
persurfaces without extrinsic curvature (i.e. surfaces of 
constant CM-time). The average energy density of the 
quark fluid on such surfaces is then given by 



"^fluid u <^ u ^Suid _ / d 3 X Tfluid e 



(9) 



and similarly for (p) . We also average over several initial 
field configurations picked according to JSJ. 

The initial condition with /i n < pe evolves smoothly 
through a crossover. For the other intitial conditions, the 
system enters the region corresponding to phase coexis- 
tence in the equilibrium phase diagram and so undergoes 
a first order phase transition. The explicit treatment of 
the dynamics of the chiral fields (in the classical approx- 
imation) in space-time allows for non-equilibrium effects 
and formation of inhomogeneities. 

Next, we determine the RMS fluctuation of the fluid 
density, Ap, induced by the propagation of the Gaussian 
initial field fluctuations JSJ through the phase transition. 
First, we determine the underlying smooth density profile 
on each t = const, time slice by averaging over the surface 
of a sphere with thickness Ar ~ 1 fm, 

(p){t,r) = Af^ 1 { d 3 x B(r + Ar-\x\) Q{\x\-r) p{t,x) , 



(10) 

with Af = 47r(Ar r 2 + Ar 2 r + Ar 3 /3). This profile is de- 
termined for each initial field configuration individually. 
While averaging (p)(t,r) over "events" (i.e. initial field 
configurations), too, would lead to seemingly larger Ae 
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FIG. 1: Evolution of the average fluid energy and baryon 
density through a crossover, and a weak and strong first or- 
der transition, respectively, for initial conditions I (top) and 
II (bottom). The densities are measured in units of nuclear 
matter saturation density po ~ 0.16 fm -3 , eo = itlnPo ~ 
0.15 GeV/fm 3 , with mjv « 0.922 GeV the mass of a nucleon 
bound in infinite matter. The fat dots indicate time inter- 
vals of « 1.5 fm/c. The phase coexistence region is shaded in 
grey; the dark area depicts the thermodynamically forbidden 
region where e < p,s p — P- 



and Ap, we are interested here in density perturbations 
on scales of order 1 fm within individual events. 

On each time slice, we then define Ap as the RMS 
deviation from this coarse-grained density profile, 



Ap 2 (t) 



Jd 3 x[p(t,x)-{p)(t,r)] 2 -(e)(t,r) 
Jd 3 x (e)(t,r) . 



(11) 



We have chosen (e)(t, r), which is obtained in a similar 
way as (p)(t, r) defined in (|10fl . as a weight in the integral 
to put more emphasis on the dense regions. Weighting 
with (p)(t, r) instead leads to qualitatively similar results. 

The time evolution of the baryon density inhomo- 
geneities for initial condition set (I) is shown in fig. [21 
We start with fluctuations of the order parameter field 
only, so that initially Ae = Ap = 0; this is to show 
the minimal degree of inhomogeneity induced just by the 
transition to the symmetry broken state. As the evolu- 
tion progresses, the fluctuations of the order parameter 
field rapidly lead to density inhomogeneities in the quark 
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FIG. 2: RMS fluctuations of the baryon density with ini- 
tial condition set (I) for crossover (narrow-dots, short-dashes) , 
weak (wide-dots, long-dashes) and strong (dash-dots) first or- 
der transitions as a function of the average baryon density. 
The fat dots indicate time intervals of ~ 1.5 fm/c. 



fluid. 

One observes that the baryon density inhomogeneities 
are sensitive to the dynamical evolution. The energy den- 
sity inhomogeneities in this model are smaller and show 
a weaker dependence on the initial baryon density and 
are therefore not shown here. 

For large initial baryon density the expansion proceeds 
through the region of first-order phase transitions. Here, 
the effective potential exhibits two local minima within 
the "phase coexistence" region of the equilibrium phase 
diagram (see e.g. fig. 1 in or figs. 2-4 in [l4j) and 
so in some region of space the order parameter can be 
"trapped" in the symmetric phase until reaching the spin- 
odal instability 15]. This effect is more pronounced the 
stronger the first-order phase transition, i.e. the smaller 
the entropy per baryon. Consequently, density pertur- 
bations can only wash out after the double-minimum 
structure of the effective potential has disappeared and 
the order parameter "rolls down" to its new vacuum. 
There is therefore reasonable hope that these inhomo- 
geneities created during the non-equilibrium phase tran- 
sition are present in the final state, contrary to those from 
the initial state. However, even for a crossover substan- 
tial inhomogeneities could be present in the final state if 
they "freeze" shortly after passing the point where V e g is 
flattest (or where the chiral susceptibility peaks, respec- 
tively) . 

Fig |3 shows our results for the set (II) of initial condi- 
tions, corresponding to fixed initial baryon density p cq = 
1.7;0o but different initial energy density e eq . Again, the 
amplitude of the density contrast is substantialy larger 
for a strong first order transition (e oq = 1.4eo) than for 
a crossover (e eq = 2.9eo). 
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FIG. 3: RMS fluctuations of the baryon density with ini- 
tial condition set (II) for crossover (dots), weak (dashes) and 
strong (dash-dots) first order transition as a function of the 
average baryon density. The fat dots indicate time intervals 
of ~ 1.5 fm/c. 



IV. DISCUSSION 

We have shown that the non-equilibrium dynamics of 
the order parameter field in heavy ion collisions can lead 
to large density inhomogeneities on the order of Ap/ po ~ 
0.1 — 1. Further, that the amplitude of the fluctuations 
depends on the structure of the effective potential: the 
effect is stronger for a first-order phase transition than 
for a crossover. 

What kind of experimental signatures could arise? By 
analogy to inhomogeneous Big Bang nucleosynthesis |23 , 
which indeed is sensitive to fluctuations of the baryon to 
photon ratio, one might expect that the relative hadron 
abundances in heavy ion collisions are modified, too. 
This is because the densities of various hadron species 
depend non-linearly on the energy- and baryon density 
of the hadron fluid and so fluctuations do not average 
out. 

The present model is too schematic to allow for quan- 
titative predictions of particle production. Nevertheless, 
experimental data for relative hadron multiplicities could 
be analyzed, for example, within the following simple 
model for an inhomogeneous decoupling surface to test 
for the presence of (energy-) density inhomogeneities. 
Thermal model fits to measured particle abundances are 
commonly performed within the grand canonical ensem- 
ble, where the density of any hadron species i can be 
expressed in terms of the temperature T and the baryon- 
chemical potential p B . Usually, a uniform temperature 
and baryon-chemical potential is assumed |23j |. On the 
other hand, to test for inhomogeneities, T and pb could 



be taken as Gaussian random variables. The average 
density of species i € {tt, K, N, ...} is then given by 



p-i(T,p B ,AT,Ap B ) = / dT P(T; T, AT) 



IX. 

x J dp B P{vb]~P B i^b) Pi(T,p B ) , (12) 



with pi(T,ps) the actual "local" density of species i on 
the decoupling surface and 



P(x; x, Ax) ~ exp 



(x — x)" 
2 Ax 2 



(13) 



the distribution of temperatures and chemical poten- 
tials. The essential point is that ~pi (T, ~p B , AT, Ap B ) ^ 
Pi(T,~p B ) if AT, Ap B =/= 0. The main contribution to the 
integrals in (fT^fl is not from T and p B but from the sta- 
tionary point of the integrand. For rare and heavy par- 
ticles, where quantum-statistical and relativistic effects 
can be neglected, there is an exponential enhancement 
of the density with (Ap B /T) 2 and (AT/T) 2 Mthis can be 
shown by a saddle-point integration of (|12fl 24] . 

From our results presented in the previous section we 
expect that AT, Ap B should be significantly larger than 
zero if decoupling occurs close to the first-order phase 
transition boundary. They should be smaller, perhaps 
nearly zero, when the dynamical trajectory did not cross 
the phase transition line. Hadron abundances at RH1C 
and SPS energies could be studied within such a model 
to search for the presence of inhomogeneities and to 
analyze their energy dependence (25|. In fact, central 
Au + Au collisions at top AGS energ y pr oduce relatively 
cool but very baryon-dense matter |2q | and could also 
probe the phase transition line j^]. Other observables 
should also exhibit some sensitivity to inhomogeneities, 
e.g. Hanbury-Brown-Twiss correlations for pions or 
production cross sections for light (anti-) nuclei, formed 
by coalescence of (anti-) nucleons [28j |. 
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